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Abstract 

This paper addresses the properties of Continuous Interior Penalty (CIP) finite element solutions 
for the Helmholtz equation. The h- version of the CIP finite element method with piecewise linear ap- 
proximation is applied to a one-dimensional model problem. We first show discrete well posedness and 
convergence results, using the imaginary part of the stabilization operator, for the complex Helmholtz 
equation. Then we consider a method with real valued penalty parameter and prove an error estimate 
of the discrete solution in the i/ 1 -norm, as the sum of best approximation plus a pollution term that 
is the order of the phase difference. It is proved that the pollution can be eliminated by selecting the 
penalty parameter appropriately. As a result of this analysis, thorough and rigorous understanding 
of the error behavior throughout the range of convergence is gained. Numerical results are presented 
that show sharpness of the error estimates and highlight some phenomena of the discrete solution 
behavior. 

Key words. Helmholtz equation, high wave number, pollution, continuous interior 
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1 Introduction 

The numerical solution of Helmholtz equation using the finite element method (FEM) in 
the medium to high wave number remains a challenge due to the strong pollution effects 
that are present in this regime. It is known that when the standard Galerkin method is 
used a so called scale resolution condition must be satisfied (see [17] ) in order to achieve 
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a quasi optimality estimate that is robust in the wave number k. Invertibility of the 
linear system also holds only under certain conditions on the relation between k and the 
discretization parameters h and p. This in particular imposes the use of high order finite 
elements and seems to exclude the possibility of using the simplest choice of piecewise 
afline elements. In this latter case the standard Galerkin finite element method has to 
be modified in order to obtain an efficient method. Such modifications often takes the 
form of least squares terms giving additional control of certain residual quantities, either 
in the element or on element faces. For low order finite elements there are a number of 
works on stabilized methods, typically using Galerkin least squares approaches and some 
results on the effect of the stabilization on the dispersion error exist in the one dimensional 
case, see [H], or for an early example of the use of face based residuals see [18]. Another 
possibility is to use discontinuous Galerkin methods and in this framework it has been 
proven by Feng and Wu [13] that provided a penalty on the jumps of derivatives over 
element faces is added to the formulation the linear system is always invertible. Similar 
results were obtained using the continuous interior penalty method in a recent work by 
Wu [22] and numerical investigations showed that the pollution error could be greatly 
reduced by choosing the stabilization parameter appropriately. For wave-number-explicit 
error analyses of other methods including spectral methods and discontinuous Petrov- 
Galerkin methods, we refer to [20l 124] . 

In the present work we continue the investigations initiated in [22], this time focusing on 
the one dimensional case and the effect of the penalty operator on the errors in amplitude 
and phase. Throughout the paper, C is used to denote a generic positive constant which 
is independent of k, h, f. C may have different values in different occurrences. We also 
use the shorthand notation A < B and B < A for the inequality A < CB and B < CA. 
A ~ B is for the statement A < B and B < A. First we will give alternative proofs 
of some of the results given in [22], showing for methods using a stabilization parameter 
with non-zero imaginary part the linear system is always well posed and the following 
error estimate holds 

\\{u-u h )'\\ ^(A^ + mina,^ 2 ))!!/!!, 

where || • || denotes the L 2 -norm. Then we consider the case when the stabilization 
parameter is real and by constructing the discrete Green's function we derive an error 
estimate where the error is written as the sum of the best approximation error and a term 
proportional to the phase error. We prove a relation between the phase error and the 
stabilization parameter and show that for a particular range of values for the stabilization 
parameter, under a mild condition on the computational mesh, the pollution error is 
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eliminated, leading to the optimal error estimate 

\\(u - u h )'\\ < kh\\f\\. 

These results are finally verified computationally in several numerical examples. 

This paper is organized as follows. In Section 2 we study the one-dimensional model 
problem and introduce the CIP-FEM. Pre-asymptotic error estimates in H l - and L 2 - 
norms are derived in Section 3 for any k > 0, h > and imaginary penalty parameters. 
In Section 4, we consider the dispersion analysis of the CIP method and obtain the phase 
error estimates between the wave number k of the continuous problem and some discrete 
wave number k^ for different real penalty parameters. The discrete global system was 
solved explicitly in Section 5 via the theory of fundamental system, it plays a major 
part in the stability and pre-asymptotic error analysis. In Section 6, the stability and 
error estimates are proved directly and we can choose appropriate penalty parameter 
to eliminate the pollution effect in this section. Extensive numerical tests are given 
in Section 7 to show some phenomena of the discrete solution behavior and verify the 
theoretical findings, and we come to the conclusion in Section 8. 

2 The model problem and its discretization 
2.1 The Boundary Value Problem 

Let Q = (0, 1) and let on Cl the boundary value problem (BVP) Lu = —f on be given: 

(1) u"(x) + k 2 u(x) = -f(x), iGll 

(2) «(0) = 0, 

(3) w'(l) - ifcu(l) = 0, 

where, for simplicity, f(x) G L 2 (Q) and k is known as the wave number. We assume that 
k 3> 1 since we are considering high-frequency problems. 

Notation 

By L 2 {VL) := H°(Q), we denote the space of all square-integrable complex- valued func- 
tions equipped with the inner product 

(v,w) := / v(x)w(x) dx and the norm := \J (w, w). 
Jq 
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We use the notation H S (Q) the Sobolev spaces of (integer) order s in the usual sense. Let 
||-|| s and |-| s denote the usual full norm and seminorm on H s , respectively. 

Existence and Uniqueness in H 2 (0, 1) 

The BVP (DO)-© has a unique solution in the space H 2 (0, 1). For the proof see, e.g., 
[3]. The existence of the solution is concluded from the following construction. 

Inverse Operator 

The Green's function of the BVP ffl-© is 



1 f sin /cxe ifcs , < x < s, 

(4) G(x,s) = -{ ' " " ' 

k [ sinks e lkx , s < x < 1. 

The solution u(x) of (JTJ) (J3j) exists for all k > and can be written as 

(5) u{x) = I G{x, s)f(s) ds, 

Jo 

and we have, 

f 1 \ coskxe iks , < x < s, 

(6) u (x) — / H(x,s)f(s)ds where H(x, s) = < 

Jo [isinkse , s<x<l. 

Lemma 2.1. The BVP ©-© has a unique solution in H 2 (0, 1) and for f E L 2 (0, 1) 

(7) N| <k-' n/ii , 

(8) Hi < 11/11, 

(9) H 2 <(! + *) 11/11 • 

Proof. See Douglas et al. [11]. □ 



Remark 2.1. T/ie aforementioned results are valid also for the adjoint problem ([[]), ([2]) 
and u'(l) + ifcw(l) = 0. 

2.2 The Continuous Interior Penalty method 

Let A4h be a uniform mesh on Cl that consists of n sub-intervals Kj = Xj), 1 < j < 

n, where Xj = j/n. Note that Xj, 1 < j < n — 1 are interior nodes and x is the Dirichlet 
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boundary node. The stepsize is h = 1/n. For the ease of presentation, we assume that k 
is constant on O. 

For any function v, denote by Vj = v(xj+) and vj = v(xj-) if the one-sided limits 
exist. We also define the jump \v\. of v at a node Xj as 

[v}j ■= vj -v+, l<j < n - 1. 

Now we define the "energy" space V and the sesquilinear form a^(-, •) on V x V as follows: 

V:={veH 1 (Q)Av(0) = 0}n J] H 2 (K 3 ), j = 1, % ■ ■ ■ , n, 

Kj£M h 

(10) a h (u,v) := (u',v') - k 2 {u,v) - \ku{l)v{l) + J(u,v) \/u,v eV, 
where 

n-l 

(11) J(u,v) := ^lh[u'] 3 [v'} 3 + 7 /i(u'(l) -iA;n(l))(t; / (l) - iAra(l)) 

i=i 

and 7 := 7 Re + Z7 lm is a complex number. 

Remark 2.2. (^aj The terms in J {u^v) are so-called penalty terms. The penalty parameter 
in J(u, v) is 7. 

^ Penalizing the jumps of normal derivatives was used early by Douglas and Dupont 
flU/ for second order PDEs and by Babuska and Zldmal JE/ for fourth order PDEs in the 
context of C° finite element methods, by Baker for fourth order PDEs and by Arnold 
^ for second order parabolic PDEs in the context of IPDG methods. More recently it has 
been proposed and analysed for fourth order PDEs by Hughes et al fTE/ and for singularly 
perturbed elliptic or parabolic problems by Burman and co-workers ^ \3$. 

(c) Notice that we here add a least squares penalty on the boundary condition as well. 
This enhances the continuity of the bilinear form and appears to be necessary for the a 
priori error estimate proposed below. 

It is clear that J(u, v) = if u G H 2 (Q) is the solution of (EE])-© and v G V. Therefore, 

(12) a h (u,v) = (f,v), VveV. 
Let Vh be the linear finite element space, that is, 

Vh ■= {vh £ H 1 ^) : Vh(0) = 0,Vh\Kj is a linear polynomial, j = 1, • • • ,n} . 
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Then our CIP-FEMs are defined as follows: Find Uh £ Vh such that 

(13) a h (u h , v h) = (f, Vh) Wv h eV h . 

We remark that if the parameter 7 = 0, then the above CIP-FEM becomes the standard 
FEM. 

The following semi-norm on the space V is useful for the subsequent analysis: 

n— 1 i 

(14) lM| lifc := (lKII 2 + £l7|/>H/)". 

3 A priori error estimate for the model problem 

In this section we will use techniques similar to those developed in [7] to derive an a priori 
error estimate that holds without any conditions on the mesh parameter and the wave 
number. We present the analysis in the one dimensional case, but the extension to higher 
dimensions is straightforward. The key observations are 

1. if the complex component of the stabilization coefficient is strictly negative (or pos- 
itive depending on the sign of the boundary condition), the formulation is coercive 
on the stabilization; 

2. if the L 2 -projection is used for interpolation in the analysis, the zeroth order term 
vanishes and the bilinear form a/ l (-, ■) has enhanced continuity properties. 

These two observations lead to an a priori error estimate on the stabilization operator 
that is optimal in h. An energy norm approach combined with a duality argument is then 
used to derive an a priori error estimate of the error in the energy norm. To simplify the 
notation in this section we assume that 7 := i7 Im the extension to non-zero real part is 
straightforward. 

Let iih : L 2 {VL) (->■ Vh be the standard L 2 -projection on Vh- It is straightforward to 
show that 

(15) \\u - ix h u\\ + h\\V(u - ir h u)\\ < h 2 \u\ 2 
and 

n 1 

(16) \J(u - n h u,u- 7r h u)\^ < |7|5(l + kh)h\u\ 2 , (V 1 |(it - 7r /l n)(x i )| 2 j 2 < h\u\ 2 . 

3=1 
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In the following we will assume that kh < 1 and neglect high order contributions in kh in 
the above approximation estimates. We first prove the continuity of a/j(-, •) on the space 
orthogonal to Vh- Let 

V ± := {v G V : (v,w h ) = 0,Ww h G V h }. 
Lemma 3.1. For all v G V 1 - and all Wh G Vh there holds 

n i 

\a h (v,w h )\ < (\J(v,v)\$ + ^l-^^/i-Xx,-)! 2 ) 2 )!^,,^)^. 

3=1 

Proof. The proof follows by observing that 

a h (v, w h ) = (v', w' h ) - \kv(l)w^{l) + J(v, w h ). 

Noting that Wh is piecewise linear and after an integration by parts in the first term in 
the right hand side we have 

n-l 

a h (v,w h ) = y^ y v(xj)[wh]j + v(l)(-ikwk~(l) + w^'(l)) + J{v,w h ). 

We conclude by applying the Cauchy-Schwarz inequality. □ 

For the stabilization operator </(•,•) we have the following stability estimate. 

Lemma 3.2. Assume that 7 Im < 0. For all Vh G Vh there holds 

\J(v h ,v h ) \ + A;|^(l)| 2 = -lm[a h (v h ,v h )} 

and for Uh solution to (TL3l) then 

\J(u h ,u h )\ + A;|n^(l)| 2 = -lm[(f,u h )]. 

Proof. Immediate by the definition of a/ l (-, ■) and (ITBl . □ 

Remark 3.1. For all^y Im < Lemma \3.2l implies existence of a unique discrete solution, 
since \J(vh,Vh)\ + k\vh(l)\ 2 is a norm on Vh- 

Combining the two previous results with the consistency of the formulation and the 
regularity estimate ([9]) immediately gives us a convergence estimate for the penalty term 
J(-, •) and the error in the right end point. 

Proposition 3.1. Let u G H 2 (Q) be the solution of (JTJ) -(EJ) and Uh G Vh be the solution 
of ( PT31) . Then there holds 

\J{u -u h ,u-u h )\ l 2 +k l i\(u- u h )(l)\ < (| 7 |l(l + kh) + |7r*)Wi||/||. 
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Proof. Let u — Uh = r\ — £h with r\ = u — n^u and ^ = — n^u. By the triangle inequality 
and the error estimate (fT6{) it is enough to consider 

Using Lemma 13.21 followed by the consistency we have 

\J(Zh,Zh)\ + kMl)\ 2 = -Im[a h (&,&)] = -Im[a A fa,&)] < M^Ca)!- 
We then apply the continuity of Lemma 13.11 to bound the right hand side, 

n i 

\j(ih,th)\ + k\ui)\ 2 < (\j(v,v)\* + i7|- i (E /i " 1 i^)iT) |J(eh ' efc)l '- 

Hence, 

n i 

(it) !■/(&, + fc*i&(i)i < i^,^ + iTr^E^ 1 !^)! 2 )"' 

then the claim follows by applying once again (fT6|) . The proof is completed. □ 

After these preliminary results we are in a position to prove the main result of this 
section. 

Theorem 3.1. (A priori error estimates) 

Let u G H 2 {VL) be the solution of (DO)-© and Uh G 14 i/ie solution of (fT3|) . luzi/i 7 Im < 0. 
Then, if h is small such that kh < 1 /or all h > and k > 1, there holds 

ll^-^n^dTl + N^mina,^ 2 )!!/!! 

and 

ll(« - V>'ll < (l7l + br 1 ) (*/* + min(l, fc 3 fc 2 )) ||/||. 

Proof. Using once again the decomposition u — Uh = 77 — by the estimate (j!5p . we only 
need to estimate the error in Consider the adjoint problem, find z G H 2 (Q) such that 

(18) -k 2 {w,z) -ikw(l)z(l) = {w,^ h ) Vto G V 
and its finite element equivalent, find 2/, G 14 such that 

(19) 2 ft ) = (w/,, tt ft - n h u) Vw h G Vh. 
By Lemma 13.21 and Proposition 13.11 z^ exists and satisfies 

\J(z h ,z h )\* = \J(z - z hl z - z h )\* < (|7|3 + |7r^)^||^!l- 
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Using the consistency of the formulation and the continuity of Lemma [37T1 followed by the 
([IE]) we get 

Uh\\ 2 = ah{£h, Zh) = a-hiv, z h) 

n i 

< (\j{v,v)\* + i7|-*(E /i " 1 W a: i)iT) |J(jKfc ' Zk)|i 



< 



(M + M^X^l/llll&ll- 

Therefore, 

(20) \\k£ h \\ < (| 7 | + M -1 )^ 2 

Next we show that ^ (| 7 | + |7| _1 )II/II- m i^ct, it follows from the definition of the 

sesquilinear form a h (-,-) that 

\\k£ h \\ 2 = -Re[a h (t h ,t h )} + {? h ,Q 

< \a h (rj^ h )\ + (J(&,&) + fcl^Cl)! 2 )^^!^^)- 1 + (kh)^)\\H h l 

where we have used an integration by parts in the second term in the right hand, i.e., 

(£Ufc) = Z:;ri 1 [a]i6^(^) + (a(l)-i^(l))6:(l)+i^i^(l)i 2 5 to derive the last inequality. 
From the continuity of Lemma 13.11 and ([16]) , ( ITT)) we conclude that 

< (1^,^)1* + N-^E/*- 1 !^)! 3 )*)!^^)!* 
^(N + M^Wll/ir. 

Therefore, 

ll^n 2 < (M + ItI- 1 ) w 2 ii/ii 2 + + fcie fc (i)i a )(| 7 |- 1 (*/i)- a + (kh)- 1 ) 

< (l + | 7 | + | 7 r 2 )(l + (A;^llfll 2 



which together with (1201) proves the first claim. 

By the definition of a^(-, ■) once again and Galerkin orthogonality there holds 

ll^ll 2 = Re[a fc (&,&)] + ||A&|| 2 < \a h (v,th)\ + ((| 7 | + ItI -1 ) min(l, k 3 h 2 )\\f\\) 2 . 
< (| 7 | + ItI" 1 ) W a ||/|| a + ((l7l + iTl- 1 ) min(l, ^ 2 )||/||) 2 . 

That is, the second claim holds. This completes the proof of the theorem. □ 
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Remark 3.2. (a) Note that the above estimate does not impose any constraints on the 
choice of the mesh size h compared to k. Both estimates exhibit the standard pollution 
term, but nevertheless the errors are upper bounded by data, independently of h and k. 
This shows that the imaginary part of the stabilization gives control of the amplitude of 
the wave. 

(b) If the penalty term on the boundary condition is removed, i.e., if J(u,v) in ( fTTl) 

n-l 

is replaced by J(u,v) := 7M M 'L'[^'h' then Theorem \3.1\ still holds. This can be proved 

3=1 

by following the analysis given in J2M/ . We omit the details. As we shall see in the next 
section, the real part of the stabilization allows us to control the phase error provided the 
stabilization parameter is chosen appropriately. 

4 Dispersion analysis 

In this section we will consider the case where 7 is a real number. Using a dispersion 
analysis we will derive precise bounds on the error in the numerical wavenumber. These 
bounds are then used to prove that a particular choice of the penalty parameter allows 
to eliminate the pollution in the one dimensional case. 

4.1 Global FE-equations and discrete fundamental system 

Let <f>2, ■ ■ ■ , 0n-i, 0n} be the nodal basis functions for the space Vh satisfying <f)j(xi) = 
5ji, the Kronecker delta, for j = 1,2, ••■ ,n and I = 0, 1, • ■ ■ ,n. Then the CIP-FEM 
solution can be spanned as: 

n 

u h {x) = ^ u h,j<Pj with u Kj = u h (xj), j = 1, 2, ■ ■ ■ , n. 
3=1 

Let — 0j, i — 1, ■ • ■ , n in (TIB"]) , the CIP formulation can be rewritten as the following 
linear system: 

(21) L h U = hF, 

where 

L h = h(a h {(j) j ,(j) i )) ,U=(u h A ,F= ((/,<&)) 

V / nxn V /nxl V /nxl 

Denote by t = kh, R = -1 - 47 - t 2 /6, S = 1 + 37 - t 2 /3, we have 
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(22) 



'-7 


R 


7 


R 


2S 


R 


7 


R 


2S 



\ 



7 



7 
/? 

/? 

7 



V 



7 

25 
/? 

7 



2S-7 



7 

i? + 2 7 



R + 2-f S-2'y-itJ 



Remark 4.1. TTie product t = kh is a measure of the number of elements per wavelength 
(of the exact solution). In particular, if the stepwidth is such that t = j for integer I then 
exactly I elements are placed on one half-wave of the exact solution. 

4.2 Discrete wavenumber and Dispersion analysis 

Recall that k is the wave number for the BVP ([I])-® and that the functions e ±lkx play 
an important role in the solution of the BVP which satisfy the equation (CQ) with / = 0. 
The discrete wave number k^ for the CIP method is defined similarly by considering the 
vector v with Vj = e h ^ h and solving the following "interior" equations: 



(23) 7^-2 + Rvj-i + 2Svj + Rv j+1 + 7^+2 = 0, j = 3, • • • , n — 2. 
Denote by th = khh, the above equations are equivalent to the equation 

(24) 2 7 cos 2 t h - 
which has the roots 



t 2 \ t 2 
4 7 + 1 + - cost h + 2 7 + 1 - - = 0, 

O / 6 



(25) 



cost 



± 



47 



6 ^ 



(1 + 1) +4 7 ^ 



Some simple calculations show that | cos t~. 



47 

< 1 < 



cost+| if 7 > -1/4 + t 2 /48 and 



I cost h | > 1 > I costal otherwise. Without loss of generality, assume | costal < 1, and 
define := t^/h and k£ := t\jh. Noting that a large | 7 | may cause a large error (cf. 



Theorem 13. 1 P and that cos t\ can not approximate cos t well (cos t\ 



27+1 
27 



^ 1 at t = 0), 



for simplicity, in the following we will assume that —1/6 < 7 < 1/6. Physically, case (— ) 
describes a propagating wave whereas case (+) describes a decaying wave [T4] . 
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Lemma 4.1. Assume that t = kh < 1, — ^ < 7 < I, then we may show 



(i) 



cost \ 



1 + - 

L ~ 2 



(ii) J/ 7 = -1/12, then - k\ < /W; 

/. .\ T « i i/i , 6 cost — 6 + t 2 cos t + 2t 2 
Hlf 17 " 7o| < ^ ^ere 7o = 12(1 _ CQst)2 : 

Proof. Denote th = % and from ff25|) . we have 

t 2 

(26) 

and 
(27) 



then \k h — k\ < /c/z. 



1 — cost; 



< 



i + f + V( 1 + f) 2 + 4 ^ 2 



t 2 2 t 2 

(1 + -) +4 7 t 2 = 4 7 (l-cos^) + l + -. 

6 D 



Clearly, t~ = k^h e (0, §) (cf. ([26])). It follows from (E5J) and (ETJ) that 



(28) cos^-l + - 



4 7 )t 4 



2(l + f + V( 1 + S +4 7 ^)(l_f + v /(l + f) ^ + 4 7 t 2 ) 
(| + 2 7 )t 4 



2 + | + 4 7 t 2 + 2 v /(l + f) 2 + 4 7 t 2 

= (1 + 6 7 )t 4 

~ 2(6 + t 2 + 12 7 (1 - cost") + 6 7 t 2 ) 

(1 + 6 7 )t 4 _ (1 + fry) (1 + 6 7 + 12^(1 - cos t^)/t 2 ) t 6 
~ 12 12 (6 + t 2 + 12 7 (1 - cost") + 6 7 t 2 ) ' 

which together with (12 6 p implies the first inequality of (i). The second inequality and (ii) 
can be proved easily as follows: the inequality sin# > -9, V6 1 G (0, |) implies 



(29) t\t^-t\<\^-t\\t^ + t\< 
and it is easy to show that: 

f2 



O • t h~ t ■ t h+ t 

2 sin — 



sin 



|cost h — cost I 



cost 



h 



2 12 



<t 6 , 



cost 



t 2 t 4 

1-- + - 

2 24 



<t 6 , 



which implies that the second inequality of (i) and (ii) hold. 

In the following, we turn to prove the last inequality. Note that cos t^ is the function of 
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7 and t h , 7 D satisfies cost h (j ) = cost and hence t h (j ) = t. By some simple calculations, 



cost h — cost 



(l + f) +4 7 t* l + |-J(l + f)%4 7o t 



47 



47o 



(l + f) +47^ -V (l + f) +4%t 2 



(1 + f + + f ) 2 + 47t 2 )(l + f + y/(l + f ) 2 + 4 7o t 2 ) 



< 



< 



4 7 t 2 



4(7 - To) 



t 2 

1+ 6 



4 7o t 5 



;i + f) +4 7 t 2 + v(i + f) +4 7o t 



t 4 < l7-7o|t 4 



and f[25j) therefore, 



* |*h - *| < |cost h - costl < t 4 I7 - j \ < t 2 h, 



which implies that (iii) holds. This completes the proof of the lemma. 



□ 



Remark 4.2. Note that the phase difference between the exact and the linear finite ele- 
ment solutions obtained is 0(k 3 h 2 ) (cf. fl\ [75]/). While for the CIP-FEM, if the penalty 
parameter 7 is close enough to ^ the phase difference is O(kh) and, as a result, the 
CIP-FEM is pollution free (cf. Theorem \6.2\ below). Figure[J\ gives a plot of the optimal 
penalty parameter 7 G versus t for < t < 1. 



-0.0835 
-0.084 

-0.0845 
-0.085 

-0.0855 
-0.086 







0.2 0.4 0.6 0.E 



Figure 1: The optimal penalty parameter versus t = kh < 1. 
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5 The discrete Green's function 

To construct the discrete Green's function, we first find the inverse of the stiffness matrix 
Lh. Inspired by the formulation of the Green's function for the BVP (cf. @), we find 
Gh = L^ 1 of the following form: 

(30) G^J^t AmM ' 



. j f% ^iA^j f% 

h ,Ti4 = e h . 

By the definition of r]i, i — 1, 2, 3, 4 there holds the facts: 



where r\\ = e lk h h ^r] 2 = e lk h h ,r]s = e lfc ^,774 = e lk h h 



(31) 771772 = 773774 = 1, Vi + V2 = 2cost h , 773 + r] 4 = 2cosi£. 
If |t| < 1/6, by some simple calculations, we can get 

(32) |cost+-l|>3. 
Without loss of generality, assume 1 774 1 > 1 773 1 , it is clear that 

(33) 1 774 1 > 3 and \r] 3 \ < -. 



From (1301) . the solution of (1211 is represented as 

n 

(34) u Kj = hY^Ghj,m(f,<f>m), J = 1,2, •••,7i, 

m=l 

and hence the CIP-FEM solution is given by 

n 

i=i 

To represent the derivative of the CIP-FEM solution, we define a n x n matrix Hh as 

(35) Hh,j,m = Gh,j,m ~ Gh,j-l,m 1 < j < Here Ghfl.m '■= 0. 

It is clear that 

n 

(36) u' h {x) = ' ] - = ^2 H h,j,m(f, 0m), Vx G [Xj-i, Xj], j = 1, • ■ ■ ,71. 

m=l 

Throughout this section let C denote a general function that may have different ex- 
pressions at different places but is bounded (uniformly) by some constant independent of 
k, h, and the penalty parameters. We first state a simple but useful lemma without proof. 
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Lemma 5.1. Suppose < t < 1, if \b\ < 0i|a|, < o\ < 1, a, b and o~\ are independent 
of the penalty parameter. Then 



(37) 



ill 



Ct 



a — bt a 

The following lemma presents estimates for H^^ m . 
Lemma 5.2. Assume that t = kh < 1, k > 1, < < §. Then 



(3* 



h,j,m 



cos(jq)e lmt * +Ct + C V l- m , j < m, 
i sin(mt^)e ij ^ + Ct + Cr)f~ j , j > m, 



where C is a general function which is bounded by some constant independent of k, h, 
and the penalty parameters. 

Proof. The proof is divided into four steps. 

Stpe 1. Solving for A m ^ and B m ^. Gh,j, m are determined by the system of equations: 

(25* — 7)Cr/ l! l, m + RGh,2,rn + lGh,3,m = &l,m, 
RGh,l,m + 2SGh,2,m + RGh,3,m + jGh,4,m = #2,771, 

lG h 

,,n—3,m ,n—2,m 

+ {2S - ^Gn 

l,m ,n,m ^n— 1,771) 

lGh,n-2,m + (R + 2^)Gh,n-l,m + {S — 2^j — 'lt)Gh,n,m = #n,m, 
k lGh,j-2,m + RGh,j-l,m + 2SGh,j,m + RGh,j+l,m + lGh,j+2,m = 3j,mi 

where 3 < j < n — 2 in the last equality of the above system and Sj tTn , 1 < j, m < n are 
the Kronecker delta. 
Formula (|2"3"|) yields 



(39) 



(40) 



7 r/r 2 + Rvr 1 + 25 + ify + 7^ = 0. 



We first consider m = 5, • • • , n — 3. From (1301 and (HOI) , the system (1391) is reduced to 
the following system of eight equations: 



(41) 



f Eti ^(25 - 7 + ^ + 7»7?)4rM = 0, 
Eti ^ 2 (^ rl + 25 + RVi + lvf)A m , = 0, 

Ell VT~ 2 [{ivf + ife?," 1 + 2S + + W^m,*] = 0, 

Eti C" 1 E(7^" 2 + RVi 1 + ZS)A m , + (R Vl + ^)B m ,i] = 0, 
Eli C [(7C 2 + RvlMm,* + (25 + ify + 7^ 2 )S m ,J = 1, 
Eti C +1 [(7C 2 )^ m , + (i&fc 1 + 25 + Rth + ^)B m ,i] = 0, 
Eti vr 1 bvi 2 + RVi 1 + 2S- 1 + (R + 2 7 )^] B m>i = 0, 
I Eti ^ [7C 2 + (R + 27)^ _1 + 5 - 2 7 - it] B mi! = 0. 
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Plugging (140]) into the first seven equations of (j4"Tj) gives 

Eti(7^ rl + ^ + 7^)^m,i = 0, 

Etl l A m,i = 0, 

(42) Eti(^ + lVi)vr(B m ,i ~ A m>l ) = 0, 

Eti(7C 2 + Rv^)vmm, - B mti ) = 1, 
Eti7C" 1 (^-^ m ,) = o, 

Eti 7(C' -2 + ^5^ = 0. 
By R = -1 - 47 - t 2 /6, 5 = 1 + 37 - t 2 /3, the eighth equation of (|41]) yields 



(43) 



E(i-j-(i + 1)*?*" 1 + 7(1 - ^r 1 ) 2 - i*)^, = 0. 

t=l 



Then, by simplifying (|4"2"]1 and (I4U1) . a 8 x 8 system which is equivalent to the system ( T4ll 
can be obtained: 



(44) 









A 




z 




v 2 _ 




B m 




_ _ 



z= [-1/ 7 ,0,0,0] J 



where A m = [A m>1 , A m>2 , A m ^ A mA ] T , B m = [B m>1 , B m>2 , B mfl , B mA ] T , and the z-th (i 
1, 2, 3, 4) column of the matrix U m , V\, V 2 are stated as follows: 



/ n7 2 \ 



( Vi 1 +Vi\ 

1 






V 2 (:,i) 



( \ 



(Li 

\ b * J 



(45) 
(46) 



m = (r/r 1 - 2 + m)rj? 7 



h = (l - \ ~ (l + j)vr X + 7(1 - V^) 2 ~ it)r£ 



1,2,3,4. 



Next we consider m = 2, 3, 4, n — 2, n — 1, n, there will be less than 8 equations, that is, 
the linear system is underdetermined, however, we can show that the system (I44p gives 
a special solution. We only prove the case m — 2, other cases (m = 3, 4, n — 2, n — 1, n) 
can be obtained similarly, we leave the derivation to the interested reader. When m = 2, 
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(47) 



from (130]) and (140]) . the system (139}) is reduced to the following system of five equations: 

Eli Vi [(25 - 7)^ + (ity + 7ifc a )3y] = 0, 
Eli »7? [i&7r%,i + (25 + ity + 7^)^2,J = 1, 
Eli Vf [(lVi 2 )A2 ti + (Rvi 1 + 2S + R Vl + Jvf)B2,i] = 0, 
Eli V?- 1 b% 2 + Rvi 1 + 2S- 1 + (R + 2 T )77j = 0, 
Eli V? bvi 2 + (R + 27)^' + 5-27-1*] Bjy = 0. 
We remark that, although the above system is underdetermined, is uniquely deter- 

mined by ( 130]) . As a matter of fact, (14"?]) can be viewed as a system of five unknowns 
B 2t i,i = 1,2,3,4 and Ei=i Vi^2,i- As we just mentioned, a solution of ( 147|) can be ob- 
tained from (jUJ) with m = 2, because of the following facts. The last three equations of 
(T4T]) are the same as the last three equations of (j4"T|) (with m = 2). The first equation of 
(147j) can be obtained from the sum of the first equation of (1411) and the fourth equation 
of (1421) (with m — 2). Similarly, the second equation of (1471) by substracting the second 
equation of (142}) from the fifth equation of (141"]) (with m = 2). 

For m = l, the system (139]) is reduced to the system of four equations: (Vi + V 2 )B 1 = z. 
In the following, we will solve (1441) . First, assuming that the matrices used are all 
invertible, implying that their determinants are not equal zero. Then, we will get A m = 
—V~ 1 V2U~ 1 z, B m = V~ 1 ViU~ 1 z, 1 < m < n, and we can also know Bi = V~ x z, where 

V = Vi + v 2 . 

Step 2. Estimating ai and bi. In order to estimate A m and B m , we prove in this 
step the following assertions: 

(48) |di| = |a 2 | < t 2 , a 3 = a^ 2 ™, |a 4 | > 6 1 774 1 ^ , 

(49) M>|f, N<J, |6 3 | < fl^l 1 "", M<^M n , 

(50) |ai& 2 — ci 2 &i I = |t 2 (?7i — "^72 ) I < 2* 2 , |«3^4 — b^a^\ < 2t 2 |?74| _n |a 4 | . 

where 774 satisfies (133]) . 
It follows from (1261) that 



|ai| = |2cost7/ - 2\ \rft\ < t 2 , \a 2 \ = |2cost7/ - 2\ \r%\ < t 2 . 
Using the identity 7/3 = 77 J 1 and (j45|) we get 

«3 = ^3 (^3 + r/4 - 2) = r/7 2n a 4 . 
It follows from ([32]) and g5]) that 

M = \v2(Vs + ^4 - 2)| = |t# | |2cost+ - 2| > 6 \r,l\ . 
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Therefore ( l4"8j) holds. 

Next, we turn to prove fl49]). Noting that < | 7 | < 1/6, from (]46j). fl3Tj) . and ([25]) 
have 

t 2 - 



\b 2 



where 



< 



t 2 i + | + + f) 2 + 4 ^ 2 



cost \ 



i + f + v / ( 1 + f) 2 + 4 ^ 2 



sin i h — t 



(I) + (II), 



(I) 



t 2 1+j + V / (l + f) 2 + 47t 2 / 1 
3 2 V 



i+f + V(i + ?) 2 + 4 ^ 2 



(l + f) 2 + 4 7 t 2 -l-f 



< 



(II) 



i + l + v / ( 1 + f) 2 + 4 ^ 2 



l-(cos^) 2 -t 



2t 2 (l-f + V / (l + f) 2 + 4 7 ^)-2t 



I _ *i 

3 12 



4 7 



t 3 



< 



(^2(1 - f + ^(1 + f ) 2 + 4Tt=) + 2) (l + f 
we therefore arrive at 

(51) l&al < (I) + (II) < 



(l + f) 2 + 4 7 t 2 ) 
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Noting that 772 = 771 , it is clear that 



bi = H(l - I - (l + ^C 1 + 7(1 - % 



1\2 



if 



2itotf. 



Obviously, |6i| > 2t - |6 2 | > I*. From ([33 



|6 3 



%" 1 



(l + -)% 1 + 7 (l- % - 1 ) 2 -i f ) 
r#(l - I - (l + ^rfe 1 + 7 (2 cos t+ - 2)»fo 1 - it 

2 7 t 2 



X "3 + 



i+i+v / ( 1+ f) 2+4 ^ 2 



773 1 - it 



< 



1 it 

3 



g ) <ofo.l 



3 1*1 



1-n 



Similarly, 



M = Kl 



t 2 

X -3 + 



2 7 t s 



l + f + V(l + l) 2 + 4 7 t 2 



--Va 1 ~ it 



3 

<2 



This completes the proof of (H9|) . 

It remains to prove (J50]l . We derive from (pE5]) - (pE6]) , fl3B . and ([25]) that 

|ai62 - 026i| = (»7i + ffe-2)(l-^-- (l + ^r/ 1 + 7 (l-r/ 1 ) 2 -it^ 

- fa + 772 - 2) (l - j - (l + |)r/2 + 7(1 - V2 ) 2 - 

= (771 + 772 — 2) ( 7 (?7i + 772 — 2) (77! — 772) — (l + (771 — 772)) 



\t 2 (vi-V2)\ < 2t 2 



Similarly, 



0364 - 6 3 a 4 = t 2 (r] 3 - 774) = t 2 a 4 r/ 4 " 



^3 + ^4-2 



t 2 a 4 ?7 4 ™ 



1-^1 



1 + 77I - 2774 ' 



hence, again from (1331) . 

I0364 — 6304] = t 2 



a 4 ?7 4 



I + 774 



< t a 4 77 4 



M + 1 

M - 1 



< 2t 2 |r7 4 r n |a 4 | . 



This completes the proof of f[50|) . 
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Step 3. Estimating A m and B m . Since t = kh < 1 and > 1, from (1331) . we have 



(52) 



Next we estimate 



detV 



By some simple calculation, we have 



det V = [(r/ 3 + 774) - (771 + 772)] [(o 2 - ai)(&4 - 63) - (62 - 6i)(fl4 - 03)] 



(53) 

where and 6j is defined by fj45|) and (|46p respectively. We analyze and estimate each 

6 2 t 

term of det V. From fl49l). it is clear that I— < -. Hence, 



(54) 61 — 6 2 = b\(l + 6*it) where 9\ is a general function satisfying \8i\ < -. 

5 

It follows from Ol. 11381). and (1511) that 



(55) (bi — 6 2 )(a 4 — a 3 ) = 6104(1 + 6 2 t) where 8 2 is a general function and \6 2 \ < -■ 

3 

From (|48])— (|49]) and (132]) . we have 

2 2 

|(o 2 - oi)(6 4 - 6 3 )| < -t 2 |o 4 | < -t|6ia 4 |. 

3 5 

It follows from (1531) . 055 j) . and the above inequality that 

det V = 6ia 4 [(% + r/ 4 ) - (77! + rj 2 )} (1 + 3 t), 
11 



where #3 is a general function and #3 < — . Therefore from Lemma 15.11 

15 

( 56 ) TTT7 = T ' 

det V 6i04<7 

where a := 773 + 774 — (7/! + 772). Note from (J2"5]l and (j3T|) that 



(57) 



7 



(i + f) w 



7(1 + Ct 2 ). 



In order to estimate Bi, we consider the first column of V*, the adjugate of V. From 
(152|) and (j4"8|) - (j5"0]) . by some calculations, we have 

/ a 3 6 4 - 6 3 a 4 + 6 2 a 4 - a 2 6 4 + a 2 6 3 - 6 2 a 3 \ / bia^Ct \ 



V*(:,l) 



—bid4 + 0164 + 6304 — 0364 + 61CZ3 — 0163 
&1CI4 — ai6 4 + a 2 b A — b 2 a 4 + a\b 2 — a 2 b\ 
\ a 3 b 2 - a 2 b 3 + a x b z - 6i<2 3 + a 2 b x - a x b 2 ) 



-6104(1 + Ct) 
6ia 4 (l + Ct) 
\ ri4 n b ia4 Ct ) 
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hence, from ([56]) and ([5 



(5 



8) B 1 = V~ 1 z = — 1 —V*z = -^—( --)V*(:,1) 

' dety det V V 7/ v 1 



( Ct \ 
1 + Ct 
-1 + Ct 

V vi n ct ) 



We turn to estimate A m and B m for m > 1. It follows from the definitions of £/ m and 
z that, 



'/l 



(59) tf- 1 * 







2—m 




\ 






-Vl)(V4 ~Vl)(V2 


-Vi) 








2—m 

V2 ViV3V* 






1 


(^3 


- mXm - m){m 


-m) 




7 




2—m 

V3 V1V2V4 










— 773) (772 -773) fat 


-m) 








9 — 777 








\ (Vi 




-Vi) 


J 



[1 + cr 



Vi - V2 

74 
% - ??4 

\ Vl~V3 / 



where we have used (l3T|) and ( j57j) to derive the last equality. 

Next we estimate V*V X . Clearly, 14(:,2) = VL(:,1), Vi(:,4) = VL(:,3), and so is V*Vi. 
It follows from pg])-(pO|) and ([52]) that, 



(60) 



V 



V*V 1 (:,[1,3]) = V*V 1 (:,[2,4}) 

( a 2 b i - a 2 b 3 - b 2 a A + b 2 a 3 0364 - b 3 a A \ 

a x b 3 — ai&4 + bid/L — bia^ 6304 — 0364 

a 2 b\ — a\b 2 a 2 bi — a\b^ — a^b 2 + a$bi 

a\b 2 — b\a 2 a 3 b 2 — a 3 bi — b 3 a 2 + b 3 a t ) 

ct_ vl n ct \ 

1 + Ct rj^Ct 

Vl n {Vi-V2)Ct 1 + Ct 
\vl n {Vi-V2)Ct - v -^(l- m Ct) J 
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From ([56]), (J59D, QSDJ, (J33j|, and = |r/ 2 | = 1, we have 



(61) 



detV 

( v? - vT 



VVxU^z 



(1 + Ct) 



Ct + Ct 



V2 - Vi 
^ ^(i + a) + Ct 



Ct \ 

sint h 

sin(mt^) + Ct 
V r/^Ct + r7 4 - 2n - 1+m C J 



Similarly, again from (E2J), (J56J) and (H5j>— <j5Dj). 

(62) (y*y 2 )([i,3],:) = -(y*y 2 )([2,4],o 



Ct 



Ct 



ru n Ct 



ru n Ct 



Vl n (Vi-V2)Ct vI n (Vi-V2)Ct V ^ n ( V4 Ct-l) -1 + Ct 
It follows from (p|h ([55]) and (|62]1 that, 

(63) A m ([l,3]) = -A m ([2,4]) : 



det^ 

„m 
72 



(T/V 2 )([l,3],:)f/-^ 



+ 



771 -^ 2 sin 

y + J 

Step 4. Finishing up. It is time to consider H h; j m . Let wj = [771,772,773,774], — 
[(771 - l)??? -1 , (772 - l)r]t\ (773 - l)7?i -1 , (774 - 1)7/I _1 ] for j > 1. From (J3Q), and 
fl63|) . we have, for m — 1, 



#M,i = G Ml i = = e u >* + Cr/J 1 + Ct = C = isin(^)e^ + + C, 

#hj,i = Gh,j,i - G h j-i tl = wjB 1 = ism(t~)e ijt ^ + Ct + CV74 
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and for m 


> 1, 




Hh,l,m 


— Gh,l,m — 


w^A m = cos(t h ' 


Hh,j,m 




Gh,j-l,m = w jA, 




Gh,j,m 


Gh,j-l,m — B\ 


Hh,m,m 




4 

Gh,m— l,m ^ ^ 
i=l 



^ '/4 > 

Jmt- + g t + Crf 4 - m , 



1 < j < m, 



4 

= W T m B m + ^(5 m>l - A^Jr}™- 1 = W T m B m 
i=l 

= isin(mt-)e im ^ + C?t + C, 

where we have used J2t=i(B m ,i — ^m^vT -1 = ® ( c ^- the sixth equation in (J12]) ). This 
completes the proof of the lemma. □ 

From Lemma 15.21 and ( |36|) , we have 

n j 

(64) u'Jz) = ^ H h,j,m(fi <t>m) = isin ( mt h) eijt *(f> 0m) 

m=l m=l 

71 71 

+ £ cos(j^)e im ^(/,0 m )+t^C ; (/,0 m ) 

m=j'+l m=l 

n 

+ Y m~ lm ~ jl C (/, <p m ), Vx e fo-i, 1 < j < n. 

771=1 

Comparing with the continuous case (EJ) we see that the first two contributions in the right 
hand side of (I64p consists of the discrete travelling wave, whereas the last two perturbed 
terms will be shown to be of the same order as the interpolation error. 

6 Stability and Pre-asymptotic error estimates for the CIP- 
FEM 

In this section, we consider the stability and error estimates of the CIP-FEM solution in 
the discrete semi-norm for real penalty parameters. 

Theorem 6.1. Under the conditions of Lemma \5.2\ the CIP method ( fl3|) attains a unique 
solution Uh that satisfies the stability estimate 

(65) K|U < ll/H . 
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Proof. Let us estimate each term in the definition of ||-|| lh (cf. (fl~4l) ). First, from ( 164)) . it 
is clear that 



\u'Jx 



< 53 I (/,^m) I < 11/11, Vx G [xj-i, Xj], j = l,>--,n, 



h 

m=l 



and hence, 
(66) 

Secondly, 
which impies 



Therefore, 



II < 



IKL'I = W h {Xj + ) -u' h (xj- )| < \u' h (xj + )\ + \u' h (xj- )| 



< 



n-1 



Em^ikl-i 2 < 

3=1 



1 

n— 1 2 



3=1 

This completes the proof of Theorem 16.11 □ 

Remark 6.1. This stability estimate for the CIP-FEM (as well as FEM) is of the same 
order as that of the continuous problem (cf. Note that the estimate holds for real 

penalty parameters in [— |, |] under the condition kh < 1 in current one-dimensional 
setting. The same result has been proved for the one-dimensional FEM in fTh^ . For 
stability estimates of the CIP-FEM for higher- dimensional problems, we refer to |HJ/ 
which, particularly, gives estimates for imaginary penalty parameters under the condition 
k 3 h 2 < 1. 

Theorem 6.2. Under the conditions of Lemma \5.2l 

(67) ||u - u fc || 1>fc < (kh + \k~ - k\) ll/H < (kh + k 3 h 2 ) 
If, furthermore, 7 = — then 

(68) \\u-u h \\ hh <(kh + k b h 4 ) 
If, furthermore, (7 - 7 Q | < then 

( 69 ) \\ u ~ u h\\ lh < kh\ 



Here 7 G is defined in Lemma 4jJ_ 
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Proof. Suppose n < k 2 , that is, k 2 h > 1, otherwise, (|69|) is proved by using the Schatz 
argument [19]. To estimate the last perturbed term in (16%|) . define go to be the largest 
integer less than or equal to — In t/ In 3. From (133|) . it is clear that 



(70) 
Define 



1 774 1 q < 3 c/ < t for g > g and go < In < k. 



^0 - 



0<x<xi, 
0, x > X\. 



Denote by Xj = for j < and Xj = 1 for j > n. We make use of the formulation of u'(x) 
in ([6]) and the characterization of u' h (x) in (1641) to obtain: For x & Kj, j = 1, 2, • • • , n, 



\u [x) - u h [x 



< 



m=0 

^ /*%rn-\-l 



„1 n 

/ H(x, s)f(s) Y] (f) m (s) ds - u[ 

j 

H(x,s)f(s)<j> (s) ds\ + Yl 

m=l 

(H(x,s)-cos(jt-)e im ^ 
\rj 4 \-^ l/l ds 



x 



%m— 1 



m=j+l Jx m-\ 



Xm + 1 



(H(x, s) - i sm(mt h )e ijt » ) f<p r , 
ds + t\ 



ds 



Xm+1 



< 

no 



/ \f\<ls + J2 \(ismkse ikx -ism(mt^)e ij ^)f(j), 

J m=l 

+ / |/|ds+ ^ / (cos/b;e ifes - cos(j^)e im ^ 

■'zj-a m=i+l "^-i 



ds 



ds 



+t 11/11 + 



J2 



l/l da 



m=l 



+ /l5 



L 2 ([xo,xi]U[a;j_2,a; : ,- + i]) 



L2 ([^i^j 2 ]) 



where ji = max{j — g — 1, 0}, ji — min{j + g + 1, n} and we have used the Lagrange 
Mean Value Theorem to derive the last inequality. Noting that (m + j) \t7 — 1\ = (m + 
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j)h \ ^h ~ ^\ —^\k h — k , the above inequality yields 

W{x) - u' h {x)\ <{t + \k- h - k\) ii/n + hh \\!\\^ M) + ^ \\f\\ LH[Xj _ 2 , Xj+l]) 
+ (go/0* ll/ILaa^^D > Wx G j = i, • ■ ■ ,n. 

As direct consequences of the above inequality, we have 
(71) ||(« - u h )'\\ 2 L2{n) <{t+ \k~ h - k\ ) 2 H/ll 2 + g 2 /* 2 ||/f + h | 

<(*+|^-fc|) ? 



where we have used q$h < t (cf. (1701) ) and h < t 2 (since k 2 h > 1) to derive the last 
inequality. 

\[(u - u h )%\ = \{u'{xj+) - u' h {xj+)) - (u'(xj-) - u' h (xj-))\ 
< \u'{xj+) - u' h (xj+)\ + \u\xj-) - v! h {Xj-)\ 

<(t+ \K - k\) ii/n + h* \\f\\ LH[X0 , Xl]) + hi \\f\\ L2{[Xj _ 3>Xj+2]) 

+ M^II/ll i2([x( ._ 1)1)a:0 . +1)2]) . 

Since I7I < 1/6, 

n-l 

(72) £ | 7 |fc |[(« - u h )'} 3 \ 2 <(t+ \k~ - k\) 2 ll/H 2 + q 2 h 2 ll/f + h ll/H 2 

3=1 

<(t+|^-A;|) 2 ||/|| 2 , 

which together with (I7ip implies fl67|) . By using Lemma [4.11 we can complete the proof 
of the theorem. □ 

Remark 6.2. (a) This theorem shows that the pollution error in H x -norm is controlled 
by the phase difference k — k^ . Ihlenburg and Babuska [15), ITb^ obtained the same result 
for the FEM in the one dimensional case. Recently, the authors 12^ showed for the 
CIP-FEM and FEM in higher dimensions that the pollution errors in H x -norm are of the 
same order as the phase difference obtained by dispersion analyses. 

(b) The pollution effect of the CIP-FEM in one dimension can be eliminated by chosen 
appropriately the penalty parameters (cf. (I69I) ). It is well-known that, the pollution effect 
exists in the FEM while in one dimension, it can be eliminated by a suitable modification 
of the discrete system but using the same stencil (cf. J^j). Note that the stencil of the 
CIP-FEM (j 7^ 0) is different from that of the FEM. We refer to [2^/ for similar results 
on discontinuous Petrov-Galerkin methods. 
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7 Numerical Evaluation 

Throughout this section, we consider the BVP with constant right hand side f(x) = — 1. 



7.1 The discrete wavenumber 



Unlike the best approximation, the CIP-FEM solution is, in general, not in phase with 
the exact solution. Numerical tests show that the discrete solution has a phase delay with 
respect to the exact solution when — | < 7 < 7 D and has a phase lead with respect to the 
exact solution when j Q < 7 < | which is similar to the FEM solution [TS]. Hence we can 
choose an appropriate value of the stabilization parameter to eliminate the phase error. 
"Optimal" values of 7 are those in a neighbourhood of 7 D . This is shown in Figure |2l 
where the real and the imaginary parts of both solutions are plotted for k — 10, kh — 1. 
There is no phase error for the CIP-FEM solution. 



it = 10, is = 10 



fe = 10, n = 10 





Figure 2: No phase error of the CIP-FEM solution with 7 = 70 for k = 10, n = 10. 

On a uniform mesh, the numerical dispersion relation of CIP method is 



(73) cost h (7) = y — , 

where t = kh. For fixed 7, the right-hand is a function of t, and is used for computation 
of the discrete wavenumber that governs the periodicity of the CIP-FEM solution. In 
Figure El the functions y\ = cost = cost^(7 D ), yz = cos tJ L {— 1/12), y 3 = cost^(O) and 
1 2/4 1 = 1 are plotted. At t c = -y/487 + 12, the functions y^ (i = 2,3) reach absolute value 
1; the numerical solution switches from the propagating case to the decaying case. The 
value t c corresponds to a cutoff frequency for the numerical solution [2T] . 
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For fixed k, the convergence k h (:= t h /h) — > k is visualized by cost h — > cost = cos kh 
as ft. — > 0. The curves begin to deviate significantly at about kh = t = 1. 



Solution of the Helmholtz Equation 

2 I 1 1 1 .- 

1.5 ■ 




'0 2 4 6 8 10 



Figure 3: Convergence of discrete to exact wavenumber via comparison of cost h (7) for 7 = j , — 1/12, to cost. 
The cutoff frequency t c = v§ f° r 7 = —1/12, t c = y/l~2 for 7 = 0. 



7.2 Error of the best approximation and CIP-FEM solution 

Consider in Figure [4] log-log-plots of the relative error := \u — Ui\i/\u\\ of the best 
approximation and the relative error e c := \u — Uh\i/\u\\ by choosing 7 = 7o for different 
k. Note that the errors first stay at 100% on coarse mesh, then start to decrease at a 
certain meshsize, and then decrease with constant slope of —1 (in log-log scale). This 
illustrates that the CIP-FEM solution is convergent to the best approximation and there 
is no pollution error for the solution. We are interested in the critical number of DOF 
where the relative error begins to decrease (see for instance [15]). We can see from Figure 
H] that the critical numbers of DOF for both the best approximation and the CIP-FEM 
solution with 7 = 7 D are about N = [-]. 

For general 7, the critical number of DOF N c can be predicted using the methods of 

\kZ - k < - w 1. 

1-3 

If 7 does not depend on t, N c follows from the Taylor expansion equation 



n c = J — '—^k') (7^ — ), N c =[ — y (7- 

V 24 / wr 12 ; V720/ v 12' 

The formula of the critical number of DOF for CIP-FEM solution is similar to FEM 
solution when 7 7^ —1/12, we consider the 7 = —1/12 case in Figure [5) It shows that the 
predicted critical number of DOF is very good, especially for large k. 
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Figure 4: The relative error of the best approximation and CIP-FEM solution with 7 = 7 in _ff 1 -seminorm and 
predicted critical numbers of DOF for k = 10, k = 40, k = 100 and k = 1000. 



Figure 5: The relative error of the CIP-FEM solution with 7 = —1/12 in iif 1 -seminorm and predicted critical 
numbers of DOF for k = 40, 100, 400. 

Figure El illustrates the relative error of the CIP-FEM solution for general 7 other 
than 7 and —1/12, say, 7 = —0.08 and 7 = — O.li, for k from 1 to 1000 on meshes 
determined by k 3 h 2 = 1. It is shown that the relative error can be controlled. For small 
k (1 < k < 50), the relative error decreases rapidly with k, for large k (k > 100), the 
relative error is dominated by the term k 3 h 2 . It verifies the estimates given by ( 1671) in 
Theorem 16.21 and Theorem 13.11 The pollution effect does exist for the two choices of 7. 

In Figure [7J the ratio e c je^ a computed with the restriction kh = 1, is plotted for k from 
1 to 1000. Obviously, the ratio (in the left of Figure [7J is increasing with k on the line. 
We remark that the ratio line in the right of Figure [7J is increasing with k and converges 
to a constant. This is due to that the relative error of the CIP-FEM solution with 7 (a 
pure imaginary number with negative imaginary part) is bounded at any range by the 
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Figure 6: The relative error of the CIP-FEM solution with 7 = —0.08 (left) and 7 = — O.li (right) in _ff 1 -seminorm 
with constraint k 3 h 2 = 1 for k from 1 to 1000. 

magnitudes of min{l, k 3 h 2 } and kh (cf. Theorem 13. ip . For large k (k > 100), the ratio 
ec/efea ^ 1 + min{l, k 3 h 2 }/kh = 1 + min{l,fc} (for kh = 1), i.e., the ratio e c /eb a < C. 
This shows that the imaginary part of the stabilization gives control of the amplitude of 
the wave. 




200 400 600 800 1000 200 400 600 800 1000 

Wavenumber k Wavenumber k 



Figure 7: The relative error ratio e c /e(, a of the CIP-FEM solution with 7 = —0.08 (left) and 7 = —O.li (right) 
to the minimal error _ff 1 -seminorm with constraint kh — 1. 



7.3 Eliminate the pollution error 

From Figure [6] and Figure [7J we know that the pollution error is present for general 7, 
but Figure [8] shows that the relative error ratio is controlled by the magnitude kh when 
we choose an appropriate parameter, say 7 = j , for n = k up to 1000. The line does 
neither increase nor decrease significantly with the change of k. 
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Figure 8: The relative error ratio e c /e-ba of the CIP-FEM solution with 7 = 7 to the minimal error iJ" 1 -seminorm 
with constraint kh = 1. 

8 Conclusion 

This paper provides some work for analyzing the dispersion and error of CIP method. We 
have show the following: 

1. The CIP method guarantees existence and amplitude control for properly chosen 
sign of the imaginary part of the stabilization operator. 

2. There is numerical pollution for general 7 and the error is mainly influenced by the 
pollution term for large k. 

3. There are many possible "good" choices of parameters to eliminate the pollution 
term. Indeed, provided kh < 1 the stabilization parameter may be chosen in an 
0(h) interval of the ideal value j . 

Future work will address the questions to what extent these results can be made to carry 
over to the multidimensional case and to higher polynomial orders. 
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